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ABSTRACT 

We examine the stability of a standing shock wave within a spherical accretion flow onto a gravitating star, 
in the context of core-collapse supernova explosions. Our focus is on the effect of nuclear dissociation below 
the shock on the linear growth, and non-linear saturation, of non-radial oscillations of the shocked fluid. We 
combine two-dimensional, time-dependent hydrodynamic simulations using FLASH2.5 with a solution to the 
linear eigenvalue problem, and demonstrate the consistency of the two approaches. Previous studies of this 
'Standing Accretion Shock Instability' (SASI) have focused either on zero-energy accretion flows without 
nuclear dissociation, or made use of a detailed finite-temperature nuclear equation of state and included strong 
neutrino heating. Our main goal in this and subsequent papers is to introduce equations of state of increasing 
complexity, in order to isolate the various competing effects. In this work we employ an ideal gas equation 
of state with a constant rate of nuclear dissociation below the shock, and do not include neutrino heating. We 
find that a negative Bernoulli parameter below the shock significantly lowers the real frequency, growth rate, 
and saturation amplitude of the SASI. A decrease in the adiabatic index has similar effects. The non-linear 
development of the instability is characterized by an expansion of the shock driven by turbulent kinetic energy 
at nearly constant internal energy. Our results also provide further insight into the instability mechanism: the 
rate of growth of a particular mode is fastest when the radial advection time from the shock to the accretor 
overlaps with the period of a standing lateral sound wave. The fastest-growing mode can therefore be modified 
by nuclear dissociation. 

Subject headings: hydrodynamics — instabilities — nuclear reactions — shock waves — supernovae: general 



1. INTRODUCTION 

The explosion of a massive star involves the formation of a 
shock wave surrounding its collapsed core. The dynamics of 
this shock wave, and the settling flow beneath it, are central 
to the explosion mechanism. A critical rate of neutrino heat- 
ing below the shock \yill drive a purely spheric al explosion 
(iBeflie & Wilsonlll985HBurrows & Goshylll993h . In practice 
this threshold may be rea ched only for low-rn ass cores (pro- 
genitor masses < 12Mq; iKitaura et al.ll2006 ^. A successful 
spherical explosion does not develop in mor e massive stars 
that form iron cores ( ILiebendorfer et al]|2001 ). 

For this reason, it has been generally believed that non- 
spherical effects are crucial to the explosion mechanism, 
at least in the case of progenitor stars more massive than 
~ 12M0. Several related instabilities of this flow have been 
identified. The first to be discovered was a convective insta- 
bility below the sho ck wave t hat is driven by neutrino heating 
jHerant et al.lll992l; see also lEpsteinlll979t) . There are also 
double-diffusive instabilities associated with the transport of 
ener gy and lepton number at or below the neutrinosphere 
(IWilson&M avle 1988; Bruenn et al.1ll995h . 

The material which flows through the shock (predominantly 
iron, silicon, and oxygen) is only weakly bound to the form- 
ing neutron star, but loses roughly half of its kinetic en- 
ergy to nuclear dissociation. The Bernoulli parameter of the 
shocked material with this dissociation energy subtracted be- 
comes substantially negative. Much of this dissociation en- 
ergy can be restored if protons and neutrons recombine into 
alpha particles, but this generally requires an outward expan- 
sion of the shock from the radius of ^ 100-150 km at which 
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it typically stalls. The successful expansion of the shock does 
not, however, require imparting positive energy to the entire 
post-shock fluid: buoyancy forces will drive a global finite- 
amplitude instability in the presence of large-sc ale density 
inhomogeneities. The stability analysis given in iThompsonI 
( 120001) shows that a non-spherical breakout of the shock is 
then possible. 

Recent two-dimensional (2D) studies of the core collapse 
problem have revealed a persistent, dipolar oscillation of the 
standing shock as an important degree of freedom in the col- 
lapse problem. In some circumstances, it has been shown 
that low-order spherical harmonic perturbations (especially 
£ = I and 2) are linearly unstable while the flow below the 
shock is still laminar A strong instability was uncovered in 
an axially symmetric hydrodynamic simulation of a standing 
shock, in the case where the fluid has an ideal, polytropic 
equation of state andjhe flow upstream of the shock is hy- 
personic dBlondin et al.ll2003l: iBlondin & Mezzaca ppa 2006). 
This corresponds to the case where the settling flow below the 
sh ock has essentially ze ro energy. A linear stability analysis 
bv lFoglizzo et al.l (1200 7'). which assumed the same ideal equa- 
tion of state, revealed high frequency modes of this instability 
to involve a feedback between ingoing entropy and vortex per- 
turbations and an outgoing sound wave which further perturbs 
the shock. We provide further evidence for this interpretation 
in our paper, particularly for fundamental, low-order modes. 

Interest in the shock stability problem arises, in good part, 
from the prominence of the dipolar oscillation that is observed 
along the axis of s ymmetry in the mos t c omplete 2D core col- 
lapse simulations (iBuras et ani2006b llal iBurrows et al.l 120071: 
iMarek & Jankall2007h . The association of such an oscillation 
with a linearly unstable mode (the Standing Accretion Shock 
Instability or SASI) should, however, be treated with caution 
when both nuclear dissociation and neutrino heating effects 
are present. A dipolar oscillation can be excited by convec- 
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tive motions (i Fernandez & Thompsonll2b08l hereafter Paper 
II), and convection generally is present below the shock near 
the critical heating rate for an explosion in the case of more 
massive progenitors. 

Given the number of competing effects that enter into the 
full core collapse problem, we have designed a set of hydro- 
dynamic simulations of a standing accretion shock which em- 
ploy equations of state of increasing complexity. In this paper, 
we use an equation of state in which the fluid flow is nearly 
ideal, with a fixed specific energy removed from the flow im- 
mediately below the shock (so as to mimic the most important 
effect of nuclear dissociation). The cooling rate is a strong 
increasing function of density and becomes significant only 
at some finite distance below the shock, which is freely ad- 
justable. In the calculations reported in this paper, we do not 
include neutrino heating, with the aim of separating the shock 
instability from the forcing effect of convection. 

We simulate the time e volution of the pertu rbed shock us- 
ing the FLASH2.5 code (iFryxell et al.ll2000h . In the hnear 
regime, we demonstrate detailed agreement between hydro- 
dynamic simulat ions and the solution to the differential eigen- 
value system of iFoglizzo et al.l (l2007l) . The linear growth rate 
is reduced dramatically as the ratio of the dissociation energy 
to the gravitational binding energy at the shock is increased, 
but we still find growing modes with several radial overtones. 
A similar trend is observed when the adiabatic index is de- 
creased. We also find that the nonlinear evolution of the shock 
is characterized by an expansion driven by turbulent kinetic 
energy at nearly constant internal energy, with the kinetic en- 
ergy reaching < 10% of the time-averaged internal energy. 
The r.m.s. displacement of the shock from its equilibrium 
position is also significantly reduced with increasing dissoci- 
ation energy. Regarding the linear instability mechanism, we 
find that the most unstable, low-frequency modes always lie 
in a region where the period of standing lateral sound waves 
overlaps with the duration of the "advective-acoustic" cycle. 

The structure of this paper is the following. In §2, we de- 
scribe our initial conditions and numer ical setup. In §3.1, we 
show that the linear stability analysis of lFoglizzo et al.l (l2007h 
yields the correct eigenfrequencies, by matching its predic- 
tions with values measured from time-dependent simulations 
to high precision. In §3.2 we use the linear stability to explore 
the effects of nuclear dissociation and adiabatic index in the 
linear phase. In §4, we explore the saturated state of the S ASI 
with time-dependent hydrodynamic simulations, focusing on 
global energetics and time averaged properties. We address 
the linear instability mechanism in §5. Our work is summa- 
rized in §6. 

2. NUMERICAL SETUP FOR TIME-DEPENDENT SIMULATIONS 

We perform one- and two-dimensional, time-dependent hy - 
drodynamic simulations with FLASH2.5 (Frvxell et al. 2000). 
This is a second-order Godunov-type Adaptive Mesh Refine- 
ment (AMR) code which implements the Piecewise Parabolic 
Method (PPM) of IColella & WoodwardI (fT984). The code al- 
lows for source terms in the equations of momentum and en- 
ergy conservation, which represent an external gravitational 
field and spatially distributed heating and cooling, as well as 
nuclear burning. 

We employ a simp lified model with no h e ating, sim- 
ilar to that used b y iBlondin & Mezzacappal (l2006l) and 
iFoglizzo et alj (|2007|) . with an ideal gas equation of state of 
adiabatic index 7, the gravity of a point mass M located at the 
origin, and a constant mass accretion rate M. We differ from 



these authors in the inclusion of an energy sink that represents 
the (partial) dissociation of nuclei downstream of the shock. 

In what follows, we describe our initial conditions and is- 
sues associated with a time-dependent calculation. 

2.1. Initial Conditions 

The initial state of our simulations is a steady and spher- 
ically symmetric flow, with a net rate of mass transfer M = 
47rr^p|v,.| from the outer boundary to the center. A standing 
shock is present at r = rso, and the flow outside this shock has 
a vanishing Bernoulli parameter 
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which ensures a vanishing energy flux through the shock. In 
the above, v^, p and p are the radial velocity, pressure, and 
density of the fluid, M = 1.3 Mi 3M0 is the gravitating mass, 
and G is Newton's constant. 

The flow upstream of the shock has a finite pressure 
pi and Mach number Mi = \vr\/csi = |v^|(7/5i/pi)"'/^, 
which are related to the upstream density by pi = 
(pi/[77W[])(M/[47rpir^()])^. There is no cooling or nu- 
clear dissociation above the shock, hence the conditions M = 
constant, b = Q, and the equation of state completely determine 
the initial flow for r > rso. 

The initial flow downstream of the shock is obtained by 
solving the time-independent continuity, Euler, and energy 
equations in spherical symmetry. The effect of cooling is 
included as a source term in the energy equation. The up- 
stream and downstream solutions are connected through the 
Rankine-Hugoniot shock jump conditions that conserve the 
flux of mass, ener gy, an d momentum across the shock (e.g.. 



iLandau & Lifshitall987h . We assume that nuclear dissocia- 
tion is a nearly instantaneous process below the shock. When 
150 km, the gravitational binding energy per nucleon at 
the shock is 12Mi .3(150km/rso) MeV, higher th an the bind- 
ing en ergy per nucleon of ^^Fe, ^ 8.8 MeV (e.g.. lAudi et alj 
l2003h . In our simplified model, the effect of dissociation is 
parameterized by a constant dissociation energy per unit mass 
e. Specific internal energy equal to e is removed from the fluid 
upon crossing the shock. To account for this process during 
initialization, the equation of energy conservation across the 
shock at r = Tjo is modified to read 
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Here the subscripts 1 and 2 label the upstream and down- 
stream flovv vari ables. This yields a compression factor 
(lThompsonl2000l) 



Pi 



(l-Al72)^ + (72-l) 
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which reduces to k ^ (7+ l)/(7- 1) for A^i 00 and e = Q. 
Increasing e increases the compression factor and decreases 
the post-shock Mach number. For example, the complete dis- 
sociation of iron into nucleons costs an energy 
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where v| = IGM/r^o is the free-fall speed. 

The rate of release of internal energy per unit volume has 
the basic form 

^o=A/?"/-". (5) 

To make contact with previous calcu lations, we have adopted 
the same parameterization as used by lBlondin & Mezzacappal 
(|2()06), a = 1.5, /3 = 2.5. This choice of exponents is con- 
sistent with cooling by the capture of non-degenerate e""" and 
e~ on free nucleons when the pressure is dominated by rel- 
ativistic particles, « pT^ cx pp^^^). The rate of cool- 
ing due to the capture of degenerate, relativistic electrons 
on protons is similarly a strong function of the pressure, 

(X Yepp^g ^ PP^^^^ where pe oc nj^ oc p'/'* is the electron 
chemical potential (the last relation assumes that the pressure 
is dominated by degenerate electrons). For simplicity we use 
only a single cooling function in this work. The constant nor- 
malization factor A determines the radius at which the flow 
stagnates, and thus the ratio r»/rso. In assessing the effects 
of nuclear dissociation on the stability of the flow, we keep 
this ratio constant and vary the normalization of the cooling 
function. 

Most of the energy loss is concentrated in a thin layer just 
outside the bound of the accretor at r = . For the adopted 
scaling of with p and p, the ratio of cooling time to flow 
time decreases with decreasing entropy. The post-shock fluid 
then undergoes runaway cooling in time-dependent simula- 
tions, and the shock collapses within a few sound crossing 
times. We have therefore implemented a cutoff in entropy in 
the cooling rate. 
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^ = ^oexp[-(V.w)'], 



(6) 



where is the cooling function without cutoff (eq. ||5|), 
i = (7- 1)"' \n{p/ p^) an entropy function, and Smm the value 
of s ?A r = r» that is obtained by evolving the flow using the 
cooling function ^ = The result is that the accreted 
fluid accumulates in the first few computational cells out- 
side r,, with a minimal modification in the outer post-shock 
flow structure. The steady stat e flow so obtained is simila r in 
structure to that calculated by Houck & ChevalieJ (1 19921) for 
hyper-Eddington accretion onto a neutron star. For our choice 
of cooling exponents, and for fixed r*/rso and 7, the main ef- 
fects of increasing e are an increase in the compression factor, 
a smaller postshock Mach number, and a steepening of the 
density profile (see Fig.fTll. In going from e = to e = 0.25v'ff, 
the normalization of the cooling function is decreased by a 
factor of 127 (keeping r,/rso and 7 = 4/3 fixed). 

We normalize the radius to the initial shock radius 
Tso, all velocities to the free-fall velocity at this ra- 
dius, Vff(rs()) = (2GM/rs())'''^, the time to the free-fall time 
'"so/vff('"so), and the density to the initial upstream den- 
sity p\. Numerical values appropriate for the stalled 
shock phase of a core collapse are ~ 150 km, 
4.8 X IO'X'3 ('"so/lSO km)-i/2 km s"'. 



Vff(rso) 



3.1Mi3''^(rso/150 km)3/2 ms, and 



Pi 



fff — 
4.4 X 



lO'MojMj 3''^ (rjo/150 km) g cm ^ (assuming a strong 
shock). Here the mass accretion rate has been normalized to 
M = 0.3Mo.3Mq s-1. 

2.2. Time-Dependent Evolution 

Our computational domain employs spherical polar coor- 
dinates, and the grid is uniformly spaced in both r and 0. 
Two dimensional calculations make use of the full range of 
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Fig. 1. — Sample initial density profiles normalized to the postshock den- 
sity p2 (upper panel), and normalized radial velocity gradient (lower panel), 
for dissociation energies in the range £/v| = [0,0.25]. (Higher curves repre- 
sent larger e.) Increasing the dissociation energy steepens the density profile 
as well as strengthening the shock compression. The velocity gradient also 
becomes stronger just above the cooling layer (white circles denote the radius 
of maximum sound speed). In this particular sequence, the normalization of 
the cooling function A is smaller by a factor of 127 when e = 0.25vjj than 
when e = 0. 

polar angles, 9 = [0,7r]. We employ a baseline resolution 
Arbase = '"so/320 and A^base = 7r/192. To better resolve the 
cooling layer, one extra level of mesh refinement is added to 
all blocks satisfying r < r*+0.1(rso-r,). To avoid excessively 
large gradients in this region and enhanced cooling due to dis- 
creteness effects, we adjust the normalization A of the cooling 
function to satisfy 



:0.995 f^-e|M, 



(7) 



where V, is the volume of the i-th cell, J^o is evaluated at the 
lower end of the cell, and the sum is performed over the post- 
shock domain. The numerical factor on the right-hand side 
is obtained empirically, and depends on the radial resolution. 
This results in an initial Mach number ^ lO"-'- 10"^ at the 
inner boundary. The default FLASH2.5 Riemann solver is 
used, which we find can support high incident Mach num- 
bers Ail < 10^. We have not witnessed the appearance of 
the odd-even decoupling instability at the shock (Ouirk 1994), 
which allows us to avoid using a hybrid Riemann solver (We 
find that the hybrid solver in FLASH2.5 has problems for 
A4i > 10 in our setup.) 

In all simulations reported in this paper, we use a reflecting 
inner boundary in the radial direction. The flow at the outer 
radial boundary (situated at 2-4rso) is given by the upstream 
steady state solution. The angular boundaries at = and 
6* = TT are also reflecting in 2D simulations. 
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Our choice of a static and reflecting inner boundary is made 
for computational simplicity. In a real core collapse, the ra- 
dius of the neutrinosphere decreases with time. The effect of 
such a dynamic inner boundary on t he linear grow th of shock 
perturbations has been examined by lScheck et al.l ('2008) in a 
semi-realistic collapse simulation. They found that it facil- 
itated growth by increasing the strength of th e velocity gra- 
dient above the neutrinosphere. In addition, iBurrows et alj 
(|2006) reported evidence for a transient S ASI instability in the 
first 100 ms of thei r full 2D co llapse simulations. The analy- 
sis of the SASI by lScheck et al. (2008), like ours, used a re- 
flecting inner boundary, albeit one positioned inside the neu- 
trinosphere. Our results provide strong evidence that linear 
growth depends on the structure of the flow between the shock 
and the cooling layer, which suggests (but does not prove) that 
a more realistic mass distribution below the cooling layer will 
not have a large impact on the growth of linear pe rturbations 
above the cooling layer ^ Yamasaki & Yamadal ( l2007i) find that 
a zero-gradient inner boundary condition (dSvr/dr = 0) has a 
quantitative but not a quahtative effect on the form of the lin- 
ear eigenmodes. 

At our baseline resolution, the accretion flow in 2D remains 
spherically symmetric for several tens or even hundreds of dy- 
namical times at the shock. To excite specific modes of the 
flow, we introduce an overdense shell in the flow upstream 
of the shock, with an angular dependence given by a Leg- 
endre polynomial of a single order £. In the absence of this 
perturbation, the flow still experiences a small startup error 
that is composed of two spherically symmetric transients: an 
outgoing sound wave due to the finite (albeit small) Mach 
number at the inner reflecting boundary, and an ingoing en- 
trop y wave due to t he initial discontinuity at the shock (see, 
e.g., lLevequdll998h . Since these initial transients do not af- 
fect non-spherical modes, we made no attempt to suppress 
them by in creasing the numerical dissipation (as is done by 
iBlondin et a l. 2003 and Blondin & Mezzacappa 2006), and 
instead use PPM in its defa u lt FLASH2.5 configuration [see 
IColella & Woodwardl(IT984l) : lFrvxeU et al] (12000) for details]. 
The spherical transients can be minimized (but never com- 
pletely eliminated) by locally increasing the resolution at the 
shock, and by simultaneously decreasing the innermost Mach 
number and increasing resolution at the inner boundary. 

To account for nuclear dissociation at the shock, we use the 
fuel+ash nuclear burning module in FLASH. In the unmodi- 
fied code, fluid labeled as fuel is completely transformed into 
ash whenever certain bounds on density and temperature are 
exceeded. A fixed energy e„ac per unit mass is released and 
added to the internal energy of the fluid in an operator split 
way. In our implementation, the fuel and ash components 
both have the same ideal gas equation of state, and are sep- 
arated by a negative energy decrement e^uc = We aim to 
choose a threshold for "burning" (that is, nuclear dissociation) 
which maintains a composition of nearly 100% fuel upstream 
of the shock, and 100% ash downstream. The composition of 
the fluid is divided this way in the initial condition, with all 
fluid injected at later times through the outer radial boundary 
being entirely /mcZ. Details of our implementation of nuclear 
burning and modifications to the original FLASH2.5 module 
are discussed in AppendixlAl 

Potential numerical instabilities result from a combination 
of a high cooling rate and a low flow Mach number in the 
cooling layer We find that the time integration is stable when 
the cutoff (|6]l is imposed on because the cooling time 
is never smaller than the Courant time. (FLASH2.5 auto- 



matically restricts the timestep so as to satisfy the Courant- 
Friedrichs-Levy condition.) Our implementation of nuclear 
dissociation also avoids introducing instabilities, since the en- 
ergy extracted from the flow in one time step is kept smaller 
than the internal energy (Appendix lAll. 

The hydrody namics module in F LASH has undergone ex- 
tensive testmg (ICalderetal.ll2002h . and so we focus our ef- 
forts on verifying our setup and on the interaction of the dif- 
ferent physics modules with the hydrodynamic solver. The 
most basic and complete test we can think of is the reproduc- 
tion of the linear growth rates of the SASI. The results are 
given in the next section. 

3. LINEAR EVOLUTION 

We now take a look at the linear growth of perturbations of 
the accretion flow below the shock. We find a strong reduc- 
tion in the growth rate when nuclear dissociation is allowed 
to take place in the post-shock flow. The time evolution of 
the perturbed flow is calculated with a direct hydrodynamical 
simula tion. We h ave also modified the linear stability analysis 
of Fogl izzo et al. (2007) to include the effect of nuclear disso- 
ciation at the shock. The real and imaginary frequencies that 
are obtained by the two methods are compared. The detailed 
agreement that is obtained for a variety of modes provides a 
stringent test of both the numerical simulation, and the solu- 
tion to the eigenvalue problem. 

3.1. Comparison of Hydrodynamical Simulation and the 
Solution to the Eigenvalue Problem 

Our method for generating a perturbation and measuring its 
growth in a hydrodynamical simulation is described in Ap- 
pendix |Bl and our method for calculating the linear eigen- 
modes of the accretion flow is summarized in Appendix ICl 
The three basic parameters of the flow are the ratio r*/rso of 
cooling radius to shock radius, the adiabatic index 7, and the 
dissociation energy e/v|-. 

It is possible to make a clean measurement of an individual 
SASI mode in a hydrodynamical simulation when the fun- 
damental is the only unstable mode (at a given £). We have 
run simulations both in regions of parameter space where this 
condition is satisfied - as predicted by linear stability analysis 
- and where it is not. The presence of unstable overtones can 
be gleaned from the time evolution of a particular coefficient 
in the Legendre expansion of the shock radius, which deviates 
from a sinusoid of exponentially increasing amplitude. 

First we reconsider the stability of the zero-energy accre- 
tion flow (e = 0), and specialize to an adiabatic index 7 = 
4/3. In the parameter range relevant to core-collapse super- 
novae, we find that the most unsta ble modes are i = \ and 
£ = 2, in ag reement with th e work of iBlondin & Mezzacappal 
(|2006|) andlFoglizzo et al.' (2007). The eigenvalue analysis of 
Foglizzo et al. (2007) correctly predicts the location of critical 
stability points of the fundamental and first radial overtone, 
for f = |0, 1,2,3|. We have also solved the differential sys- 
tem of iHouck & Chevalied (|I992|) . finding that not only the 
critical stability points but also the growth rates do not agree 
with what we measure in our si mulations (although results for 
l=Q are identical to those of .Foglizzo et al.ll2007h . In the re- 
mainder of this subsection, we refer to simulation results for 
which only the fundamental is unstable. 

The left panels of Fig. |2] show growth rates and oscilla- 
tion frequencies^ for the modes t= 1,2,3, in the case where 

^ We have chosen to plot eigenfrequencies in units of the inverse free-fall 
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Fig. 2. — Growth rates ajg,ow (upper panels) and oscillation frequencies ojosc (lower panels) of the SASI, in units of inverse free-fall time at the shock, versus 
the ratio of stellar radius to shock radius r»/rso- In the left panels, the unperturbed flow has zero energy flux (e = 0) and upstream Mach number M[ ~ 87. In 
the right panels, specific internal energy e = 0.25vjj is removed from the flow just below the shock; and Mi = 5. The lines give the solution to the eigenvalue 
problem of Foslizzo et al. (2007) for different spherical harmonics: £ = i (red), £ = 2 (green), and £ = 3 (blue), with thick lines representing the fundamental 
mode and thin lines the first radial overtone. The dashed lines show the effect of neglecting the cutoff in the cooling function (eq. (S)). (See Appendix Icl for 
details.) Symbols show the eigenfrequencies obtained from time-dependent hydrodynamic simulations: for £ = I (down-triangles), £ = 2 (squares), and £ = 3 
(up-triangles). Uncertainties in the fitted parameters are smaller than the symbol size. (See Appendix |B]for details.) 



the flow upstream of the shock has a high Mach number 
Ail 87. There is very good agreement between the cal- 
culated eigenfrequencies and the output of the hydrodynamic 
simulation: 1 -3%, 1 -2%, and 1 -2% for the real frequen- 
cies of the £ = 1, 2, and 3 modes, respectively, and 2-4%, 
2-8%, and 5 - 8% for the growth rates. Our method for es- 
timating the uncertainty in the measured mode frequencies is 
detailed in AppendixlB] The most important systematic errors 
arise from the discreteness of the mesh, the discrete summa- 
tion involved in the Legendre projection, and the discrete time 
sampling. The error bars are comparable to or smaller than the 
size of the symbols in Fig.|2] typically Sujtff ~ 5 x 10"-'. 

The agreement between the two methods of calculating the 
growth rates becomes worse at larger £. We attribute this to 
the better sampling of the lower-£ modes by the grid, which 
results in weaker numerical dissipation. (The effects of nu- 
merical dissipation in PPM typically depend on the ratio of 
cell width to wavelength: Porter & Woodward 1994.) 

When the dissociation energy is increased to e = 0.25v'ff, the 
agreement between the two calculational methods is reduced 
a bit, to 4-5% for the oscillation frequencies for both £ = I 
and £ = 2, and 18-30% and 20-42% for the growth rates. 




Fig. 3. — Growth rates (left) and oscillation frequencies (right) as a function 
of resolution, for the £ = 3 mode in a flow with r* / r^o = 0.7 (see panels a,b 
of Fig. [3. The dashed line in each panel shows the solution to the hnear sta- 
bility problem, and symbols represent measurements from hyd rodynamical 
simulations at different radial and angular resolutions. See j|2.2l for a descrip- 
tion of the grid spacing, and Appendix|B]for an explanation of the error bars. 
The dotted line shows the best fit power law to the growth rate difference for 
a radial and angular resolution increase: Ao^gmw on (ArA^)"" '^^" '". 



time tff, since this is directly related to simulation time when both the stellar 
mass and the shock radius are held fixed. 



respectively."* A more extensive exploration of the influence 



^ The absolute value of the discrepancy between the two sets of growth 
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Fig. 4. — Growth rates (left) and oscillation frequencies (right) of the fundamental £ = l (red) and £ = 2 (green) modes as a function of rt/r^o, for dissociation 
energies e/vjj = 0,0.05,0.1,0.15,0.2,0.25,0.3, and 0.33. Panels (a) and (b) show eigenfrequencies in units of the inverse free-fall timescale %, with decreasing 
curves corresponding to increasing e. Panels (c) and (d) show the same curves, but in units of the inverse of the time required to traverse the postshock cavity at 
the postshock speed, (rjo -''»)/ |v2 1. In panel (d), increasing dissociation makes the peak of the curves move to the right. Other parameters of the sequence are 
7 = 4/3 and M[ — + oo. See text for description. 



0.15 - 



0.05 - 





-0.05 



Fig. 5. — Radial overtones of the £ = I mode as a function of nuclear 
dissociation energy e. The bold line denotes the fundamental mode, and the 
thin lines the first 11 overtones. Other parameters are r^/rso = 0.4, 7 = 4/3, 
and Ml —> 00. Even though nuclear dissociation significantly decreases the 
growth rate, there is always an unstable overtone. 



rates is similar to that found for £ = 0, namely (10 ^ — 10 ^)«jj' . 



incident Mach number 
Fig. 6. — Growth rates (left) and oscillation frequencies (right) of the 
fundamental £ = I mode as a function of incident Mach number, for e / Vjj = 
0,0.025,0.05,0.075,0.1,0.125,0.15,0.175, and 0.2, with decreasing curves 
for higher dissociation energy. Other parameters are rt/r^o = 0.4 and 7 = 4/3. 
Eigenfrequencies are modified significantly when A4i < 5. The normaliza- 
tion of the cooling function A is bigger by a factor 4 — 5 for Mi =2 relative 
to A^l = 100 throughout the range in e. 

of finite e and a reduced upstream Mach number on the mode 
frequencies is made in the next subsection. For now, we em- 
phasize the good agreement between analytics and numerics. 
We have also checked the effects of resolution. Fig.[3]shows 
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Fig. 7. — Radial overtones of the £ = I mode as a function of adiabatic index 7, for a zero-energy accretion flow [e = 0, panels (a) and (c)] and a flow with 
internal energy e = 0.25\j^ removed below the shock [panels (b) and (d)]. The upper row shows results in units of the inverse free-fall time %, whereas the lower 
row shows eigenfrequencies in units of the inverse of the time required to traverse the shock cavity at the postshock velocity, (rso-''»)/|i'2l- The bold line denotes 
the fundamental mode, the thin lines the first 6 (a,c) and 1 1 (b,d) overtones. Growth rates generally decrease toward smaller 7, but there is always an unstable 
overtone. This qualitative result is independent of the size of the density jump at the shock. Other parameters are ff/rso = 0.4, e = 0, and jVfi 00. 



results for r^/rso = 0.7 and ^ = 3 [see panels (a) and (b) of 
Fig. 121. The growth rates show the expected behavior, with in- 
creasing agreement between the two methods with increasing 
radial and angular resolution. We observe a greater sensitivity 
to the angular width IS.6 of the grid cells than to their radial 
width. The real frequencies show a much weaker trend. In 
sum, a four-fold decrease in Ar and IS.9 improves the agree- 
ment in the two sets of growth rates from 7% to 0.8%, while 
the disagreement in the real frequencies stays nearly constant 
within uncertainties at 1.5%. 



3.2. Dependence of Eigenfrequencies on Nuclear 
Dissociation and Adiabatic Index 

We now turn to the effect of nuclear dissociation on the os- 
cillation frequencies and growth rates. As is observed in the 
right panels of Fig.|2] there is a substantial reduction in growth 
rate when the dissociation energy is raised to e = 0.25v|- 
(and the Mach number upstream of the shock is reduced to 
M. \ =5). One also observes a significant drop in the real fre- 
quency. This effect is further illustrated in Figs. |4^,b for the 
fundamental i=\ and i = 2 modes, as the rate of nuclear dis- 
sociation is gradually increased. 

Both trends can be partially explained in terms of an in- 
crease in the advection time below the shock. A finite nuclear 
dissociation energy at the shock increases the compression 



factor K (eq. |3]l and decreases the post-shock Mach number 



M2 



K 7H ^ 



-1-1/2 



(8) 



The advection time scales as r^Q jv^ oc k, whereas the lateral 
sound-travel time has a weak dependence on k. As a result, 
the advection time increases relative to the lateral sound travel 
time. The growth rate of a mode of fixed I peaks when the ad- 
vection time is comparable to the period of a lateral sound 
wave within the settling zone between the shock and the sur- 
face of the accretor (see ® . The peak growth rate therefore 
moves to higher r^jr^i^, for fixed I, or to a lower Legendre in- 
dex for fixed r^jr^^^, as e is increased. This can be observed in 
Figs.|2]and|4^,c. 

The advective-acoustic cycle yields a WKB growth rate 
scaling as Wgrow ~ lii|2|/fQ^ where Q is the efficiency o f the 
cycle and fg ^ 27r/a>osc its duration (Foglizzo et al.ll2007i) . For 
larger dissociation energies, fg is dominated by the advection 
time, and fg increases somewhat faster than linearly^ in k. 
The oscillation frequencies depend mainly on the flow time, 
as can be seen by scaling them to |v2|/(rso-r,) (Fig.|4}l). On 

For large values of e, the density profile steepens relative to p ~ and 
thus the velocity decreases faster than linearly with radius. The advection 
time )12t is then dominated by the innermost part of the flow. 
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the other hand. Fig. |4}; shows that the growth rates decrease 
with increasing e even after this first-order effect is removed. 

The growth rate also depends on the coefficients for the 
conversion of an outgoing sound wave to an ingoing entropy 
or vortex wave at the shock; and for the Hnear excitation of 
an outgoing sound wave by the ingoing mode. As is shown 
in Appendix|D] the first coefficient increases with increasing 
dissociation energy. In the Hmit of strong shock compression, 
it is given by eq. ( |D9t 



2(7-l)(l-A^2) 
7^2(1+7-^2) 



1)k 



1/2 



(9) 



We therefore conclude that the reduction in growth rate sig- 
nals a decreasing efficiency for the conversion of an ingoing 
entropy-vortex perturbation into a sound wave. 

To understand why this happens, note that the radial flow 
time increases relative to the lateral sound travel time as e is 
pushed upward. This makes it more difficult for the pressure 
perturbation at diametrically opposite points on the shock to 
maintain the correct relative phase. As a result, the peak in 
the growth rate is forced to larger values of r,/rso. We also 
observe that the velocity gradient at the base of the settling 
flow becomes stronger with increasing shock compression: 
the gradient scale Vr/{dvi /dr) is reduced by a factor of ^ 2 
as e is raised from to Q.25v^. We suspect that this second 
effect has a weaker influence on the growth rate, given our 
observation that that rapid growth depends on an approximate 
equality between the radial flow time and the lateral sound 
travel time in between the cooling layer and the shock (see 
®. 

The fundamental mode can be stabilized for sufficiently 
high e, as is apparent in Fig.|4] Nonetheless an unstable mode 
can always be found among the higher radial overtones. Fig- 
ure |5] shows the growth rates of the first 1 1 overtones of the 
l=\ mode as a function of e, for r^/rs,Q = 0.4. 

Our ability to find li nearly unstable modes stands in con- 
trast to the analysis of lYamasaki & Yamadal (l2007h . which 
employed a single, more reaUstic equation of state, and found 
no unstable modes in the absence of neutrino heating. But 
the comparison between both studies is made difficult by the 
fact that we are keeping the ra tio r^, /r^n fixed, whereas in the 
absence of heating Yamasaki & Yamada ( ,20 07) obtain a very 
small shock radius. It should also be emphasized that our un- 
stable modes undergo a weak non-linear development: the po- 
sition of the shock is perturbed only slightly when e > 0. 15 v| 
(see ®. Neutrino heating therefore plays a crucial role in 
maintaining large-amplitude oscillations of the shock for re- 
alistic values of the dissociation energy e (Paper II). 

A reduction in the upstream Mach number has a measur- 
able effect on the frequency and growth rate of the SASI (see 
Fig.|6]l. When A^i > 10, one reaches the asymptotic, strong- 
shock regime, and the eigenfrequencies are essentially con- 
stant. Lower values of A4 1 push up the post-shock velocity 
and the real frequency of the mode (at fixed value of r^/rso). 
Growth is substantially reduced for a weak shock when e is 
close to zero, but is only modestly affected when e is large. 

The dependence of the growth rate (£= 1) on the adiabatic 
index 7 is shown in Fig. |7]for e = and 0.25vg, in units of 
fff"' as well as |v2|/(rs()-r,). Reducing 7 has some similar 
effects to increasing e at fixed 7: the flow slows down below 
the shock, and the density profile steepens. The growth rate is 
therefore reduced, and the peak growth rate is found at higher 
radial overtones (see ®. When the radius of the shock is 



allowed to vary, the peak growth rate moves to a higher value 
of / as 7 is decreased. 

4. NONLINEAR DEVELOPMENT IN AXISYMMETRY 

The non-linear development of the SASI is addressed in this 
section. Our hydrodynamical simulations enforce axisymme- 
try, and ignore the effects of neutrino heating. Our main in- 
terest is in how nuclear dissociation limits the growth of the 
SASI. We first discuss the transition from a laminar accretion 
flow to a turbulent, quasi-steady state. We then explore some 
time-averaged properties of this state, and the global energet- 
ics of the flow below the shock. 

4.1. From Linear Instability to Saturation 

The oscillation of the shock, and the fluid below it, remains 
coherent during the initial linear phase. Consider the £ = I 
SASI mode. The shocked fluid in the expanding pole has an 
increased energy density relative to the surrounding fluid. As 
a result, the expanding pole is the source of a sound wave that 
propagates laterally to the other pole. When the radial advec- 
tion time is comparable to this lateral sound travel time, this 
pressure enhancement reaches the opposite pole just as the 
phase of the oscillation has reversed, thereby reinforcing the 
ou tward expansion of the shock. T he standing wave observed 
by iBlondin & Mezzacappal (12006') can be interpreted in this 
way (instead of the purely acoustic phenomenon proposed by 
these authors). Another visualization of the effect (for e = 
can be seen in Movie 1 in the online material. 

The non-linear development of the SASI, in the absence 
of nuclear dissociation or in the presence of neutrino heat- 
ing, involves relatively large shock deformations involving 
non-radial flows and the development of strong shock kinks 
(e.g, Blondin & Mezzacappa 2006; Scheck et al. 2008). On 
the other hand, when nuclear dissociation takes a significant 
fraction of the accretion kinetic energy and the SASI is not 
forced by convection, shock oscillations saturate at a much 
lower amplitude, with shock kinks not forming. Nevertheless, 
nonlinear coupling between different modes still takes place 
and the shock thus acts as a generator of vorticity. The transi- 
tion to the nonlinear saturated state for three different values 
of e is shown in Movie 2 in the online material. 

In the longer term, the axisymmetric flow reaches a quasi- 
steady state, with a broad range of oscillation frequencies. 
One could draw a parallel between this behavior and that 
of confined 2D turbulence, which reaches a quasi-steady 
state in which the vorticity accumulates on the larg est spatial 
scales, decaying on the very long viscous timescale (iDavidsonI 
l2004l) . The nonlinear sat uration prope r ties of the SASI have 
been studied i n 2D by Ohnishi et al.l (l2006h and in 3D by 
llwakami et al.l (|2008). Both of these studies include a num- 
ber of effects, among them neutrino heating and a particular 
semi-realistic equation of state. We defer an examination of 
the effects of heating to Paper II, and focus here on the effects 
of nuclear dissociation on this quasi-steady state. 

4.2. Time-Averaged Properties 

Since none of our simulated flows explode, we are able to 
average the properties of the fluid over a fairly long period 
of time (T ^ 10^- lO-' dynamical times at the unperturbed 
position of the shock). In a turbulent system that displays a 
quasi-steady state beha vior, a time aver age is equivalent to an 
ensemble average (e.g. lDavidsonll2064l) . The mean and rm.s. 
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fluctuation of a quantity /(r, 0, f ) are defined in the usual way, 

(10) 



(/(r,0)) = i / fir,e,t)dt, 
i Jo 



and 



Af(r,0)^{{f)-{f)') 



1/2 



(11) 

Fig. m shows the velocity field and sound speed, averaged 
over time and angle 9, in accretion flows with e = and 
0.15v|-. At any given moment one can define minimum and 
maximum shock radii, which allows the construction of time 
averages (r,_min), (r,,,nax) and fluctuations Ar^^min, Ar^^m^- 
One can distinguish four different zones in the flow, from the 
inside out: (1) the cooling layer, extending from the inner 
boundary to roughly the point of maximum of (c,); (2) the 
adiabatic envelope, which is bounded above approximately 
by (r, mill) - Arj min; (3) the shock oscillation zone, extending 
from this radius out to (rs.max) + Ars ^ax; and (4) the super- 
sonic accretion flow. It is immediately evident that, as ex- 
pected from the quasi-steady state behavior, accretion pro- 
ceeds almost the same as in the unperturbed case, with (vv) 
and (vg) closely following the unperturbed velocity profile. 
The fluctuations Avv and Avg are comparable in magnitude 
and much larger than the mean flow, although they remain 
below (cj). As expected, the angle average of (vg) vanishes. 
Notice also that the sound speed in the mean flow is every- 
where lower than in the initial configuration, a point that we 
explore below. 

A dramatic feature of Fig. [8] is the reduction in the ampli- 
tude of the shock oscillation as the effect of nuclear dissocia- 
tion is introduced into the flow. It is important to understand 
the extent to which the oscillations seen in core collapse simu- 
lations are the result of an intrinsic fluid instability, as opposed 
to a coupling of the shock to the convective motions that are 
maintained by neutrino heating. To this end, we have run a se- 
ries of simulations with e increasing from to 0.25 v|., and an- 
alyzed the change in the character of the turbulence. Since the 
growth rate of the SASI is strongly reduced by nuclear disso- 
ciation, we have chosen a ratio of accretion radius to shock ra- 
dius for which the growth rate of the fundamental peaks when 
e = (namely r»/rso = 0.4). This is also compatible with our 
choice for studies of explosion hydrodynamics when neutrino 
heating is included (Paper II). A fraction 1 -2£/v^ of the ac- 
cretion energy at the shock is available for exciting turbulent 
motions below the shock. This corresponds to 50%- 100% of 
the accretion energy for the flows that we are examining. 

The results of this study are summarized in Figs. |9] and [TOl 
and Table [U The first figure shows the amplitude of the £ = 
and £ = 1 Legendre coefficients ao,i of the shock radius, as a 
function of time, for a few runs. Their amplitude drops dra- 
matically as e increases above 0.15 v|-. Some configurations 
display significant intermittency in the oscillations, as exem- 
plified with the case e = 0.1 v|. 

The following figure shows the time-average (ao.i) and 
rm.s. Aflo.i of the Legendre coefficients ((oi) ~ 0). These 
quantities have been normalized to the time-average of the 
steady shock position for a ID run with the same parameters. 
This procedure allows us to remove a small offset from rso 
in (flo) fliat is caused by the finite initial velocity at the inner 
boundary of the simulations, as can be seen in the upper pan- 
els of Fig. |9] The £ = Q mode, being stable, settles to a steady 
value within a few oscillation cycles. For the lower values 
of e, a modest monopole shift in the shock radius remains in 
2D, as compared with ID, but this becomes insignificant for 



> 




sO 

Fig. 8. — Time-averaged profiles of velocity and sound speed in the 
fully developed, nonlinear phase of the SASI. Upper panel represents a 
zero-energy flow (e = 0), and lower panel a flow in which internal energy 
e = 0. ISvjj is removed below the shock. Notice that the shock moves over a 
smaller range of radii in the second case. Spherical and temporal averages 
are given by the solid lines, and rm.s. fluctuations by the dashed lines, for: 
radial velocity (blue), radial velocity in the unperturbed flow (lower black); 
meridional velocity (red), sound speed (green) and sound speed in the un- 
perturbed flow (upper black) Vertical lines denote: the inner boundary (left 
dotted); the inner excursion of the shock (middle dotted); the outer excur- 
sion of the shock (right dotted); the rm.s. range of the £ = component of 
the shock radius (dot-dashed). See the text for definitions of these quantities. 
Other parameters are r*/rso = 0.4, 7 = 4/3, A4i = 87. 

£ > 0.15vff. The time-averaged £ = 1 coefficient nearly van- 
ishes when half of the accretion energy is removed by disso- 
ciation (e = 0.25 v|.). 

In a real core collapse, the protoneutron star contracts by 
more than a factor of two in radius while the shock is still 
stalled. We have therefore run a series of simulations with 
different r^/rso, to check whether the envelope size has a sig- 
nificant influence on the properties of the saturated state. Fig- 
ure[TT]shows the ^ = 0, 1 Legendre coefficients for these runs, 
which have vanishing dissociation energy. As the envelope 
size increases by a factor of two, the fractional £ = Q expan- 
sion of the shock nearly doubles, increasing from 1.15(r5)iD 
to 1 .29{rs) ID, where (rj) id is the steady shock position in ID. 
On the other hand, the rms amplitude of the £ = I oscillations 
only increases by ~ 20%. 

The partitioning of the energy of the fluid below the shock 
into different components is explored in the lower panels of 
Figs. [To] and [TT] Since our ID runs do not display any grow- 
ing mode, and settle into a well-defined steady state configura- 
tion, we use them to calculate reference values of the kinetic, 
internal, and gravitational potential energies. These quantities 
are denoted by £'kin,iD, £^int,iD, and /igiav,^, respectively, with 
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200 400 600 
time [tjj,] 

Fig. 9. — Amplitudes of ^ = (upper) and £ = l (lower) components of the 
shock radius, versus time. Panels from left to right represent flows with in- 
creasing dissociation energies: e = 0, 0. Ivjj, and 0.2vjj. The horizontal lines 
bound the r.m.s. fluctuation of each Legendre coefficient. These configura- 
tions correspond to those in Movie 2 (online material). 



TABLE 1 

Time- Averaged Energy" of Flow below Shock (2D) 





e/4 


l(£grv>| 
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{Ekin.r) 


(^kiii.e) 


0.4 





62 ±4 


40±3 


1.8±0.6 


1.7 ±0.5 




0.05 


74±3 


45±2 


1.6±0.7 


1.7±0.6 




0.1 


95 ±3 


55±2 


1.4 ±0.7 


1.4 ±0.7 




0.15 


129.1 ±0.9 


74±1 


0.4 ±0.2 


0.4 ±0.2 




0.2 


200.0 ±0.7 


107.6 ±0.6 


0.2±0.1 


0.1±0.1 




0.25 


320.3 ±0.7 


159.6 ±0.3 


0.21 ±0.02 


0.13 ±0.03 
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39 ±3 


26±2 


1.2±0.5 


1.1±0.4 


0.375 





70±6 


46±4 


1.8±0.7 


1.6±0.5 


0.25 





123 ±7 


80±4 


3.1±1.0 


2.5 ±0.7 


" Ener 


gies 
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in units 


of piv2,J 


« 3.4 X 



1 O^'^Mo.sMJ ''3- (r^o / 1 50km) 1 /2 



erg. Error bars denote the r.m.s. fluctua- 



the spatial integral being carried out between the radius at 
which Cs peaks and {rs)\D- The analogous quantities (/ikin.iD), 
{Emtio), and (£'grav.2D) are obtained by integrating over vol- 
ume and averaging over time. * Table[T]displays the time av- 
eraged energies and their rms fluctuations for both sequences 
of runs. By increasing e and thus the compression rate at the 
shock, the mass in the adiabatic layer increases and so does 
the gravitational binding energy. The internal energy also in- 
creases with increasing dissociation energy, because near hy- 
drostatic equilibrium needs to be maintained. 

The internal energy remains almost constant, within fluc- 
tuations, relative to the ID time averaged case, with some 
increase for e = and e = 0.2. Its ratio to the gravitational 
energy, however, is decreased relative to the ID case by about 
10% for e = 0, about the same fraction accounted for by the to- 
tal kinetic energy. So our results point to a scenario in which 
the post-shock envelope expands as a result of the onset of 
turbulence, keeping its internal energy nearly constant. The 
fraction of the accretion energy that is converted to turbulent 
kinetic energy decreases with increasing e, from ^ 10% of 
the internal energy for e = to about 0.1% for e = 2. From Ta- 
ble[T]one can see that the total turbulent kinetic energy is very 
close to the fluctuation in the internal energy. As regards the 
sequence for different /r^o, we note that the fraction of the 

* The integration over volume is performed from the radius at which the 
angle and time average of Cj peaks to the instantaneous shock surface. 
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Fig. 10. — Time-averaged properties of the shocked flow. Horizontal axis 
labels the fraction of the gravitational energy that is available for exciting os- 
cillations of the shock. Upper panel: r.m.s. fluctuation of the £ = and £ = 1 
components of the shock radius, with points denoting mean values. Radius 
has been scaled to the equilibrium shock radius in the ID solution (which 
is stable). Middle panel: various components of the energy of the 2D flow 
below the shock, averaged over time, relative to the ID flow. Black symbols 
denote internal energy and red symbols denote ratio of internal to gravita- 
tional energy. Error bars represent the r.m.s. fluctuation. Lower panel: time 
average of the fluid kinetic energy below the shock, relative to the internal 
energy (red symbols) and gravitational energy (black symbols). Notice the 
dramatic drop in these ratios with increasing dissociation energy. The system 
expands mainly as a result of the growth of turbulent kinetic energy, with a 
smaller net change in internal energy. 

internal energy going to turbulence decreases with increasing 
envelope size, but the absolute value of the turbulent kinetic 
energy more than doubles. 

5. ON THE LINEAR INSTABILITY MECHANISM 

The nature of the linear instability mecha nism has been 
clearl y demonstrated in the WKB regime dFoglizzo et alj 
l2007h : the growth of higher radial overtones is due 
to the "advectiv e-acoustic" cycle first described by 
iFogUzzo & Tagged f2000). This cycle involves a non- 
spherical deformation of the shock, which creates entropy 
and vortex waves in the downstream flow. These perturba- 
tions become partly compressive as they are advected toward 
the star, thereby creating an outgoing sound wave that inter- 
acts with the shock. There is growth if the secondary shock 
oscillation is larger in amplitude than the initial perturbation. 
The instability is aided if the flow is strongly decelerated 
somewhere below the shock. There is no purely acoustic 
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Fig. 1 1. — Same as Fig. 1 101 but now varying Tt/r^o for fixed e = 0. Note 
that the vertical scale varies between the panels. 

instability. 

The linear stability analysis of iFoglizzo et alj (l2007h 
found approximate agreement with the real frequencies and 
growth rate s measured i n the h ydrodynamical simulation of 
iBlondin & Mezzacappal (l2006h - which involve low-order 
modes of the shock - and good agre e ment for the £ = modes 
in particular. lYamasaki & Yamadal (l2007h solved the eigen- 
value problem in an accretion flow with a realistic equation 
of state and strong neutrino heating. They also found that the 
eigenfrequencies most closely matched the advective-acoustic 
cycle. Hydrodynamical simulations show that the oscillation 
period scales with the advection time at mod erate to large 
shock radius (" Ohnishi et aT]|2006l: IScheck et alJlIO OS). 

Here we provide further insight into the mechanism driv- 
ing the fastest growing, low-frequency modes of a spherical 
shock. We show that growth involves the radial advection of 
entropy and vortex perturbations, and the lateral propagation 
of sound waves in the settling flow below the shock. This ex- 
plains the observation by Blondin & Mezzacappa (2006) that 
the communication of pressure perturbations by lateral sound 
waves plays a role in the SASI, but disagrees with their infer- 
ence that the mechanism may be purely acoustic. 

To this end, we calculate the growth rate as a function of the 
size of the settling region by solving the eigenvalue problem 
as formulated by iFoglizzo et all (I2007D for several spherical 
harmonics in a zero-energy accretion flow (e = 0). The basic 
configuration of the flow is the same as described in ^ By 



changing the ratio of to rso, we are able to change the ra- 
dial advection time relative to the time for a sound wave to 
propagate laterally within the flow. 

The upper panel of Fig. [12] shows the growth rates of the 
fundamental (thick lines) and first radial overtone (thin lines) 
of the £ = 1-4 modes, as a function of r*/rso. The lower 
panel compares the advection time of the flow from the (un- 
perturbed) shock in to r, , 

fadv= / i—r> (12) 



with the period of a meridional sound wave of Legendre index 
£. The acoustic period is plotted at two different radii: right 
below the shock, 

(13) 



1,2 



and at the base of the settling flow. 



27rr* / r» 



3/2 



(14) 



Here 2 is the sound speed downstream of the shock, and 

Cs.* = Cs,2('"so/'"*)'''^ is a good approximation to the sound 
speed at the top of the cooling layer. The self-gravity of the 
accreting material is assumed negligible, and most of the post- 
shock region is essentially adiabatic for the chosen cooling 
function. Hence, fmax.f corresponds to the longest possible 
period of a lateral sound wave of spherical harmonic i out- 
side the cooling layer. The shortest acoustic period (^ fmin.f) 
is found at the base of the settling flow. 

For any given £, the peak of the growth rate of the funda- 
mental mode is always found where tj„in,£ < fadv < tmaxj'^ the 
growth rate of the first radial overtone peaks where fmin.f < 
fadv/2 < fmax.^- The basic result does not depend on whether 
one uses the radial advection timescale, or the sum of the 
advection time and the radial sound travel time (which is 
significantly shorter). Fig. [13] focuses on the radial over- 
tones of £ = 1, which have significant growth at large val- 
ues of /"so/r,. The peak of the n-th radial overtone satisfies 

fmin,l <fadv/(n+l)<f max, 1 ■ 

The main conclusion that we draw from Figs. [T2]and[T3]is 
that the radial flow time controls the overtone of the fastest- 
growing mode, while the period of the meridional sound wave 
controls the angular order This interplay between the two 
timescales points to the advective-acoustic cycle as the mech- 
anism driving the instability. 

The £ = mode deserves special comment. We find that 
its oscillation period is nearly twice the duration of the ad- 
vective acoustic cycle, 2(fadv + fs.up), where fs.up is the time 
taken for a sound wave to travel radially from r^, to Tsq. This 
is demonstrated for a wide range of shock radii in Fig. [14] 
A purely spherical perturbation of the shock will excite a 
downgoing entropy wave. When the shock oscillation reaches 
its maximum or minimum radius, the sign of the pressure 
perturbation below the shock is generally opposite to that 
of the ent ro py p er turb ation, as m ay be seen by combining 
eqs. ini, iDl, Ids), and ||D9l. When the entropy wave 
reaches the cooling later, the sign of the change in the (neg- 
ative) cooling rate per unit volume (eq. lO) is the same as 
the sign of the entropy perturbation (at constant pressure), 
SJiffo/^Q = -[(7- l)/'j]{P-a)SS. As a result, a positive pres- 
sure perturbation at the shock^ generates an increase in the 

' As measured at its unperturbed position. 
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Fig. 12. — Upper panel: grow rates of linear modes of the shock, versus 
r* / rso for various spherical harmonics: £ = I (red), £ = 2 (green), £ = 3 (blue) 
and £ = 4 (orange). Thick curves denote the fundamental mode; thin curves 
the first radial overtone. Lower panel: various timescales in the flow below 
the shock: advection time fadv from shock to star; period fniax.f of lateral 
sound wave just below the shock; and period /niin.f of lateral sound wave just 
above the cooling layer Dotted lines show the intersection between fajy and 
'max.f I and dot-dashed lines the intersection between i^jy and fmin.f- Growth 
of the fundamental mode is concentrated where ^niin.f < fadv < 'max.£: and 
similarly < 'adv/(" + 1) < fmax.f for the radial overtone. 

cooling rate at the base of the settling flow, and therefore a 
negative pressure perturbation that is communicated back to 
the shock on the radial acoustic time. This secondary pressure 
perturbation is in phase with the pressure perturbation at the 
shock if the mode period is twice the period of the advective- 
acoustic cycle. This analysis of the £ = mode helps to ex- 
plain the importance of lateral acoustic waves in the growth 
of non-radial perturbations. 

6. SUMMARY 

We have investigated the linear stability, and non-linear 
evolution, of a standing shock in a spherical accretion flow 
onto a compact star, in the context of core-collapse su- 
pernovae. Our analysis combines two-dimensional, time- 
dependent hydrodynamic simulations with a linear stability 
analysis. The equation of state used in these calculations al- 
lows for a loss of internal energy to nuclear dissociation below 
the shock, but is otherwise that of an ideal gas. We do not in- 
clude heating effects associated with charged-current absorp- 
tion of neutrinos on protons and neutrons, or the recombina- 
tion of alpha particles. A companion paper explores how the 
threshold for an explosion is influenced by these additional 
effects (Paper 11). 

Our main results are as follows: 




Fig. 13. — Same as Fig. 1121 but now focusing on the radial overtones of 
£ = 1 at large r^/r^f,. We denote by f. nl, n2, n3, and n4 the fundamental and 
first fo ur o vertones, respectively. Other quantities have the same meaning as 
in Fig.[l2] 
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Fig. 14. — Period of the 1 = fundamental mode (black curve), which 
is very close to twice the advection time (red) and twice the duration of the 
advective-acoustic cycle (green). Parameters of the flow: 7 = 4/3 and e = 0. 

1 . — The growth rate of the Standing Accretion Shock In- 
stability (SASl) is significantly reduced when ~ 20-50% of 
the gravitational binding energy of the flow is absorbed by 
nuclear dissociation. A lowering of the adiabatic index of 
the flow (or, equivalently, a steepening of the density profile 
above p{r) ~ r~^) has a similar effect. The lower growth rates 
appear to result from a weaker coupling between the ingoing 
entropy-vortex wave and a sound wave that acts back on the 
shock. This reduced coupling is caused in part by a lengthen- 
ing of the radial flow time in comparison with the time for a 
sound wave to propagate laterally below the shock. Increasing 
the dissociation energy also pushes the strongest instability to 
lower £ at fixed r^,/rso, and to higher n at fixed i. 

2. — By solving the eigenvalue problem, we do find some un- 
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stable modes for all values of e. By contra st, the linear stabil- 
ity analysis of lYamasaki & Yamadal (l2007l) found no unstable 
modes for radial overtones < 2 at vanishing neutrino luminos- 
ity (but evidence for very strong growth of n = 3 modes, which 
we do not find), albeit with small shock radii since their cool- 
ing function was fixed. 

3. — The nonlinear development of the SASI in two dimen- 
sions is characterized by a quasi-steady state that is reached 
after the linear modes have saturated. As the dissociation en- 
ergy is increased, the fluctuations of the shock about its unper- 
turbed position are greatly reduced in amplitude. We infer that 
strong neutrino heating is needed to drive large-amplitude, 
dipolar oscillations of the shock in a realistic core collapse 
environment. 

4. — The expansion of the shock that we observe appears to 
be driven mainly (but not entirely) by the action of turbulent 
pressure; we observe a smaller absolute positive change in 
the internal energy of the flow below the shock. The turbulent 
kinetic energy saturates at < 10% of the internal energy. A 
large expansion of the shock therefore is possible only if the 
flow has nearly vanishing total energy. 

5. — The strongest growth of shock perturbations of spheri- 
cal harmonic t and radial overtone n is encountered when the 
radial flow time across a distance ^ (r-r*)/(n+ 1) equals the 
period l-nrjlcsir) of a lateral sound wave, at some radius r 



in between the cooling layer and the shock. This provides a 
compelling argument as to why the linear instability involves 
a feedback between a radially propagating entropy-vortex per- 
turbation, and a laterally propagating acoustic perturbation. 
6. — The period of a spherical (^ = 0) mode is nearly twice 
the radial flow time, that is, approximately double the period 
of the fundamental non-radial modes. This is due to a change 
in sign between the pressure perturbation at the shock, and 
the pressure perturbation that is induced in the cooling layer 
at the base of the accretion flow by the advected entropy per- 
turbation. 
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APPENDIX 

A. NUMERICAL TREATMENT OF NUCLEAR DISSOCIATION IN FLASH2.5 

Nuclear dissociation is implemented by means of the fuel+ash nuclear burning module in FLASH2.5 (iFryxell et al.ll2006 h. To 
prevent undesired dissociation effects upstream of the shock, we replace the density and temperature thresholds for burning in the 
default version of the code with a threshold in Mach number: burning takes place so long as the fluid has a Mach number lower 
than Alburn = 2. This has the effect of localizing the nuclear dissociation right behind the shock. We find a very limited amount 
of incomplete burning whenever the shock is strongly deformed, so that it runs nearly parallel to the upstream flow. Of course, a 
similar phenomenon is expected in the full collapse problem (incomplete dissociation behind a strongly oblique shock leading to 
the formation of cold and fast accretion plumes); so we have made no attempt to completely eliminate this effect. 

In principle all of the kinetic energy of accretion can be absorbed by nuclear dissociation downstream of a shock: the Riemann 
problem has a limiting solution p2 ^ oo and V2 0, with p2V2 = pivi, which for a strong shock yields e — ^ Vj/2. However, 
the numerical evolution of the problem in discrete timesteps limits the dissociation energy to be smaller in magnitude than the 
specific internal energy of the fluid, otherwise negative pressure results. In particular, since we want to maintain a steady shock 
as background flow for stability calculations, the maximum dissociation energy that can be removed in a single timestep at r = 
is the internal energy of the postshock flow after dissociation has been subtracted. That is, if we allow for the density contrast k 
to increase according to eq. Q, then we require 

e<- '^rr^a^, (Al) 

where M2 is the post-shock Mach number, eq. ([8]l. The maximum single-step dissociation energy e^-^x is obtained by equating 
the two sides of eq. ( lAll ). One finds emax — 0.213v^ for 7 = 4/3 and A^i 00. 

In reality, burning of the shocked fluid occurs within a layer of a finite width. To relax the above limit on e, we have modified 
the default fue I + ash submodule of FLASH2.5 to allow nuclear dissociation to occur in stages: instead of burning all the/Me/ into 
ash in a single step, we only allow the burning of a fraction 1 /nbum at a time. The value of «bum is adjusted empirically so as to 
avoid numerical problems at large expansions of the shock, where e may approach the local value** of v|/2. In most cases, the 
burning is spread across ^ «bura cells behind the shock. We have checked that the mode frequencies of the flow are insensitive to 
the particular choice, as long as the single-step constraint is met. 

B . PROCEDURE FOR OBTAINING EIGENFREQUENCIES FROM A HYDRODYNAMICAL SIMULATION 

Here we describe the procedure to obtain Unear eigenfrequencies from our time dependent numerical simulations, and the 
associated error. To excite a particular mode in our simulations, we add an overdense shell upstream of the shock with the 
desired angular dependence and a sinusoidal radial profile (amplitude ^ 10 %, width ^ 0.2rso). This generates an initial shock 
displacement with very small amplitude (negligible compared with the spherically symmetric displacement described in § I2.2l i. 
followed by a growth of the SASI mode on the advection timescale. 

* As we discuss in § 4, the shock oscillations saturate at a low amplitude when e > 0.15v|r, which means that this effect is negligible for the simulations 
reported in this paper. 
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To measure the growth rate Wgrow and oscillation frequency Wosc of a mode associated with a given Legendre-£, we project the 
shock surface onto the corresponding Legendre polynomial, obtaining a time-dependent Legendre coefficient, 

ae(t)= / Rs{e,t)Pe(0) sine de, (Bl) 



2 JO 

where Pi is the Legendre polynomial and Rs(0, t) is the shock surface, defined as the locus of points with pressure Pshock = \/P\P2- 
For the strong shock simulations, this has a numerical value /?shock — O.Olpivg^. 
The time-dependent Legendre coefficient is then fitted with the functional form 

at{t) = c\+C2 sin ( Wosc? + C3) + C4f, (B2) 
where Wgiow is the growth rate, Wosc the oscillation frequency, and the c,- are constants. This fitting function works well whenever 
there is only a single unstable harmonic, usually the fundamental, for a given I (see, e. g., lower-left panel in Fig. |9]during the 
first ~ lOOfff). The fitting is performed with a Levenberg-Marquardt algorithm jPress et al. 2006). 

For the fitting we choose a temporal range that begins when the Legendre coefficient starts a clean sinusoidal oscillation, and 
which ends when aiit) = (rso-r,)/10 for l>\, corresponding to a shock displacement of 10% of the unperturbed value at the 
poles. Our temporal sampling interval is fff/2, about a quarter of the radial sound crossing time from r» to Tsq. 

To estimate the uncertainty in the fitted eigenfrequencies, we assign a "measurement error" to at{t), which we compute as 
follows. The error in R^{9,t) is taken to be one-half the size of the baseline resolution, Arbase/2. The size of the angular cell 
A6'base enters through the computation of the integral in eq. (IB lb as a discrete sum. Errors adding up in quadrature then yield 

5a, = Arbase A(?base PjiOi) sin^ 6i , (B3) 

where the sum is performed over all the angular cells. This result is then used as the input error in the Levenberg-Marquardt 
algorithm. The error bars shown in Fig. |3]are the 1-sigma errors that output from the fitting routine, multiplied by 3. 

C. LINEAR STABILITY ANALYSIS WITH NUCLEAR DISSOCIATION DOWNSTREAM OF THE SHOCK 

We outline our approach to calculating linear perturbations of a spherical accretion flow with a standing shock wave. The 
surface of the star imposes a hard inner boundary to the flow, and therefore introduces a free parameter /rso, the ratio of the 
stellar radius to unperturbed shock radius. The structure of the flow is further defined by the equation of state and the cooling 
function. We generalize the ideal gas equation of state to include a finite, uniform dissociation energy e downstream of the shock, 
as is described in ^2.1J We explicitly account for the low-entropy cutoff in the cooling function (eq. ||6l). 

iFoglizzo et al.l (120071) formulate this eigenvalue problem assuming a constant-7 ideal gas equation of state. Generalizing 
this formalism to account for a constant dissociation energy e is straightforward when thermal energy is removed from the 
flow only immediately below the shock. The background flow solution is modified, but the algebraic form of the perturbation 
equations is not. We therefore use the differential system given in eqns. (10)-(13) of iFoglizzo et al.l (|2007|) . which employs the 
variables /= v,(5v, + 2q(5q/(7-1), h = S(pvr)/ pv,-, SS = (Sp/ p-jSp/ p)/(-f- 1), and SK = ■ (V x Su}) + £(£+l)(cyj)6S, 
where a; = V x v. All quantities are projected onto spherical harmonics Y". A complex eigenvalue oj = Wqsc + 'Wgiow is obtained 
by enforcing Sv,- = at r = r» + 10~'*(rso-r*). The eigenfrequencies obtained are nearly insensitive to the radius at which this 
boundary conditions is applied, so long as it lies inside the cooling layer. 

A subtlety in the treatment of the boundary conditions at the shock is worth discussing. When the shock is perturbed to a 
position + Ar and velocity Av, the equation of energy conservation across the discontinuity becomes 

1... a.a2, c?, _ 1 2, 7 f P2 , 5p2 P2Sp2\,de^ de 



-(VI - Av)^ + ^ = :7 (V2 + <5v,- Av)^ + ^iL± + ^-L±^\+ Sp2 + —Sp2. (CI) 
2 7-I 2 l-l \P2 P2 P2 pi J op dp 

Here, as before, 1 and 2 label the upstream and downstream flows in the frame of the accretor. When e = constant, the algebraic 
form of this equation is unchanged from the case e = 0. Although the Bernoulli parameter b of the background flow below the 
shock is reduced, the perturbation / = 5b does not receive additional terms. The boundar y conditions on /, h, 6S, and 6K therefore 
have the same algebraic form as eqns. (B10)-(B12), (A6) and (B15) of Foglizzo et al.l (|2007) . By contrast, the shock boundary 
conditions on the perturbation variables dp2, 5p2 and (5v,- do acquire additional terms resulting from the changing compression 
ratio K (eq. [S)). We have checked that these additional terms cancel out in the boundary conditions on /, h, 6S, and 6K. 

To account for our entropy cutoff (eq.|6]l in the linear stability calculation, one needs to add an additional term to the perturbation 
to the cooling term in the energy equation. 



5 ' 



(p-1)— + 2^ 



(C2) 



pv / \pv 

The precise analytic eigenfrequencies are somewhat sensitive Smin, so this parameter needs to be the same in both linear stability 
and simulation for proper comparison. The difference between including and excluding the entropy cutoff in the linear stability 
can be seen from Fig.|2l where results obtained without entropy cutoff are shown as dashed lines. 

D. COEFFICIENT RELATING THE AMPLITUDE OF AN INGOING ENTROPY WAVE TO AN OUTGOING SOUND WAVE AT A 

STANDING SHOCK: EFFECTS OF NUCLEAR DISSOCIATION 

We focus here on the WKB limit, where the wavelength of the excited modes is small compared with the density scale height 
below the shock. The perturbed equations of mass, momentum and energy conservation across the shock are 

pi (vi - Av) = ip2 + 6p2)iv2 + Sv, - Av), (Dl) 



ACCRETION SHOCK WITH NUCLEAR DISSOCIATION 



15 



py{vi - l^vf + Pl = {P2+ 5p2)(V2+ 5Vr- l^vf + P2+5p2, (D2) 

and eq. ( ICll ) with de/dp = de/dP = 0. Since the end result becomes fairly complicated in the case of non-vanishing dissociation 
energy, we now simplify the discussion by taking the limit Mi ^ oo. Working to linear order in the perturbation variables, we 
divide up them up into contributions from an outgoing sound wave - (moving oppositely to the background flow), an ingoing 
sound wave +, and an entropy perturbation S. Then 



P2 



5pf = c]25pt; 5vf = TCsi^—\ (D3) 



552 = -^^; (5pf = ^vf = 0. (D4) 

7-1 Pi 

Equation ( ID2b becomes 

^Pl=—, vi^Pl--, —7^^Pl = -7TTTT-:^^Pl-7TTTtv^^Pl (D5) 



where M.2 = \ v2\/cs2 =1/ V7(k-T)- Substituting for 5p2, equations (ICll l and (ID lb become 
and 

, (i-M2) , - M2 ,s 
'(1 + A^2)^" (I + M2) 



(p2-pi)Ay = 2 ^^ ^ ^ c,2Sp2- , . . , Cs2Sp2 (D7) 



respectively. One also has 

Vl=-C,2(^X2+^ ), (D8) 



for arbitrary k > 1 . Combining these equations gives 

6pl 2(7-l)(l-X2) 



(D9) 



Sp2 A^2(l + 7-^2) 

This result can be expressed in terms of the perturbed Bernoulli variable [f^ = -(7- 1)"' c^2 ^p\l Pi\ = (1 i M2)c^2^Pt / P2^ 

f 2 



/± 7W2(1+7A^2)' 
When the dissociation energy vanishes, M.2 = [(7- l)/27]'/^ and eq. ( ID9I ) reduces to 



(DIO) 



Sp?, AM2 

= —. (Dll) 

(5p2 I+27W2 

In this particular case, we find agreement with the coefficient derived in iFogUzzo et alj (l2005h (their equ ation [Fll]). That 
previo us result cannot, however, be used when £ ^ 0, because its derivation depends on the Hugoniot adiabatic (iLandau & Lifshita 
[l987h . One can also show that, in the absence of dissociation, 5p2 = 5p2+5p2 + 5p2 = 0. This is no longer the case when e ^ 0, 
because a net density perturbation arises due to the change in the compression rate [eq. ©J, 5p2 = piSti, with 

Sk = -^Av. (D12) 
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